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ABSTRACT 

Stellar black hole (BH) binaries are one of the most promising gravitational wave (GW) 
sources for GW detection by the ground-based detectors. Nuclear star clusters (NCs) located 
at the centre of galaxies are known to harbour massive black holes (MBHs) and to be bounded 
by a gravitational potential by other galactic components such as the galactic bulge. Such an 
environment of NCs provides a favourable conditions for the BH-BH binary formation by the 
gravitational radiation capture due to the high BH number density and velocity dispersion. We 
carried out detailed numerical study of the formation of BH binaries in the NCs using a series 
of N-body simulations for equal-mass cases. There is no mass segregation introduced. We 
have derived scaling relations of the binary formation rate with the velocity dispersion of the 
stellar system beyond the radius of influence and made estimates of the rate of formation of 
black hole binaries per unit comoving volume and thus expected detection rate by integrating 
the binary formation rate over galaxy population within the detection distance of the advanced 
detectors. We find that the overall formation rates for BH-BH binaries perNC is 10“^°yr“^ 
for the Milky-Way-like galaxies and weakly dependent on the mass of MBH as T oc 
We estimate the detection rate of 0.02-14 yr“^ for advanced LIGOWirgo considering several 
factors such as the dynamical evolution of NCs, the variance of the number density of stars 
and the mass range of MBH giving uncertainties. 
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1 INTRODUCTION 


The gravitational wave (GW) is the propagation of ri pples of 
the sp ace-time curvature with speed of light. Although lEinsteinl 
( Il916h predicted the existence of gravitational radiation (GR), the 
GWs have never been directly detec ted yet. From 30 years ob¬ 
servations, IWeisberg & Tavlon i2005tl found that t he bin ary pul¬ 
sar PSR 1913-1-16, discovered by Hulse & TavloJ jl974ll . exhib¬ 
ited the decrease of the orbital period and the amount of decrease 
exactly coincide with the prediction of general relativity. There 
are several types of astronomical objects that can produce signifi- 
cant amount of gravitational waves: core-co llapse supemovae (e.g., 
iMueller & Jankall997l:lYa kunin et al.l2010h . spinning neutron stars 
(e.g., lAn dersson et alj |201 ih . compact binary coalescences (e.g., 
iKalogera et alJl20M l. supermassive black ho les (SMBH) (e.g., bi¬ 
nary SMBH me rger. [B erentzen et alj 2009|) ( extreme mass ratio 
inspirals, EMRI, iHopman & AlexandeJ 2006 jMemtt_et^ 1201 ih 
and cosmological density fluctuations (e.g., Ananda et alj 2007 1. 
Among those, the compact binary coalescence (CBC) involving 
neutron stars (NS) or stellar mass BH is the primary targets for 
the second generation of ground based detectors such as advanced 
LIGO and Virgo. Up to now, only ten NS-NS binary pulsars have 
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been discovered in our Galaxy including PS R 1913-1-16, and a half 
of the m will merge within a Hubble time JO’Shaughnessv et al. I 
l2010h . Based on these binaries and theoretical studies, the rate of 
NS-NS merger within the horizon of the adva nced detectors is es - 
timated to lie between 0.4 and 100 per year iAbadie et al.ll20l5j . 
These NS-NS binaries are thought to have been evolved from pri¬ 
mordial binaries. Obviously other types such as BH-NS and BH- 
BH binaries can form by the evolution of binaries of high mass 
stars, but their rates are highly uncertain especially because of the 
lack of observed systems. Since the detection of the GWs coming 
from the CBC event will provide us with the information regarding 
the masses of the binaries, we should be able to distinguish differ¬ 
ent types of binaries through the observations of GWs in the future. 


Compact binaries can also be formed dynamically in stellar 
systems: three- body processes and di ssipative two-body processes 
by tidal effe ct dUe e & Ostriked 198 61) or GW emissi on (hereafter 
GR capture. iHansenl 19721 : Quinlan & Shapirolll987h . Stellar sys¬ 
tems such as globular clusters (GC) and nuclear star clusters (NC) 
at the galactic nuclei provide good environments for the dynami¬ 
cal formation of compact binaries. In GCs, when the core is dense 
enough, compact binaries can be formed by three-body encoun¬ 
ters. These binaries become more compact through close encoun¬ 
ters with other stars and is eventually kicked out from GCs when 
their orbital separations become very small. Some of the ejected hi- 
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narie s merge within a Hubb l e time in galacti c field dBaneriee et alj 
1201(1 IPowning et al.ll201ll : ISae et alll20l4) . On the other hand, 
binary formation by three-body processes is supp ressed by the ex¬ 
istence of massive BH (MBH) in galactic nuclei dBaumgardt et alJ 
l2004ah . Instead, coalescences of primordial compact binaries can 
be driven by the orbit couplin g with central MBH, as known as 


the Kozai effect i Kozai 


_ _ _ 196^. and could occur within Hubble 

time dAntonini & Peret^ 2012h Compact binaries, especially BH- 
BH binaries, also can be formed by GR capture in NCs due to 
the hig h stellar density and velocity dispersion at the vicinity of 
MBH dO’Learv et alj|200^ . The black holes are expected to be 
highly concentrated in the central parts through the dynamical fric¬ 
tion against low mass stars, providing a favorable environment for 
the close encounters between black holes. Such captured compact 
binaries usually have large eccentricities with small pericentre dis¬ 
tance, and thus, their waveforms will be different from the ones 
produced by the circular binaries dAbramovici et al.|ll992h . 

Numerous authors made estimates of the detection rates of 
GWs from CBCs with present- and planned- GW detectors using 
various methods: population synthesis models for primordial bi¬ 
naries dKalogera et ai.ll200 4[Belcz vnski et al]|2007li . Monte-Carlo 
simulations for GCsj[Dowmin2_etay|20H ) and Fokker-Planck sim¬ 
ulations for NCs do Leary et al.ll2009l) . and estimated that a few 
tens of events will be detect ed by new generatio n GW detectors ev¬ 
ery year (see, for example, lAbadie et al.ll201(ih . However, most of 
these studies are based on simplified models and assumptions on 
the evolution of the stellar systems and the binaries. There remain 
substantial uncertainties because of difficulties in accurately mod¬ 
eling of the systems with the large number of stars. Direct A^-body 
approach is impossible for realistic systems, and therefore, statis¬ 
tical approaches such as Fokker-Planck models and Monte-Carlo 
simulations have been used, so far. In this paper, we focus on the 
binary formation of BHs by GR capture in NCs by direct A-body 
simulations. Although, we cannot use realistic number of stars, we 
try to deduce important information regarding the binary forma¬ 
tion and evolution based on scaled-down and simplified version of 
A-body simulations. Our study should provide useful comparison 
with the statistical models. 

This paper is organized as follows. In §§2 and 3, we introduce 
the numerical method and the model for star clusters in galactic nu¬ 
clei with a central MBH. The dynamical evolution of our model is 
presented in §4. In §5, we describe binary formation in NCs and es¬ 
timate the merger rate per galaxy. In §6, the expected detection rate 
for new generation GW detectors is estimated. To give the interpre¬ 
tation for binary coalescences and their waveforms, we implement 
post-Newtonian approximations on the two-body motions. An ex¬ 
ample waveform of a binary BH coalescence in Milky-Way-like 
galaxies is provided in §7. Finally, we summarize in §8. 


2 INITIAL MODELS 

NCs are very dense stellar s ystems located at the nuclei of galaxies , 
regardless of the type (e.g.. lCarollo et alJl 19971 : lBbker£taljj200j: 
C6^et all 1 20061) . Their typical mass is l lWalcher et al.l 

2005h. the size of NCs is compar able to that of galactic GCs 
tebker et~^l2004l : ICote et alj|2006l) . thus average stellar density 
of NCs is much higher than GCs. It is als o well known that most 
of ga l axies host MBHs at the centre (e.g.. lKormendv & Richston^ 
ll995l : lFerrarese & Fordl200^ . The coexistence and correlation of 
M BHs and NCs at the centr al region of galaxies have been studied 
bv iGraham & Snitleil l l2009ll . 


Observational (e. g.. [Schbdel et'^l2009h and theoretical (e.g., 
iBahcall & Wolll Il976h studies for star clusters with the central 
MBH showed that the density and velocity dispersion increases 
steeply toward the MBH, following Kepler’s profile (i.e., cr oc 
^- 1 / 2 ) modeling of stellar systems with a central black 
hole has been done by n umerous authors (e . g., lYound [l98^: 
Goodman & Binned j^84 (Quin lan et al.l 1 19951 : ISigurdsson et al.l 


I 995 I : Hollev-Bockelmann et alj 2002fi ~ In order to generate N- 


body realizations for NCs w ith MBH, we adopt ‘adi abatic growth’ 
of the MBH as suggested bv ISigurdsson et air(ll995h . The MBH is 
assumed to grow with time as 


M'MBH(f) = 


1 ATmbh 

3 f 

\ ^MBH / “ \ ^MBH / 

(ATmbh 



t ^ Imbh 

t > Imbh 

( 1 ) 


where Mmbh and Imbh are the final mass of MBH and the black 
hole growth time scale, respectively. During the growth of the 
MBH, the stellar system is adjusted against the potential of the 
MBH. The MBH is fully grown after Imbh, and the gravitational 
potential of the MBH is assumed to follow that of the Plummer 
model 


fl^MBH 


GMmbh 
+ Embh 


( 2 ) 


where £mbh is the soft ening parameter f or avo iding unexpected 
effect at the singularity. ISigurdsson et al.l l ll995h noted that fiviBH 
should be larger than the half-mass dynamical time to ensure the 
adiabatic growth of the MBH. While they used the Hernquist model 
as the initial density distribution, we used Plummer model with 
half-mass dynamical time fdyn.i /2 ~ 2.46 for standard A-body 
scaling (i.e., G = = —4i7 = 1). 

According to recent observations (e.g.. lSch6del et alj2009li of 
the centre of the Milky Way (i.e., the vicinity of Sgr A*), the ve¬ 
locity dispersion of stars at about 1 parsec scale is nearly flat. This 
implies that the NC is almost in the isothermal state. However, it is 
not possible for the isolated stellar systems to become fully isother¬ 
mal. In order to realize nearly isothermal sphere we need a very 
large and massive system. Instead we model the NC as a isother¬ 
mal sphere surrounded by external potential which is assumed to be 
fixed. Such an external p otential could be provided by the galactic 
bulge. f^on et al.l ( 120111) have investigated a self-gravitating stellar 
systems embedded in an external potential well. They considered a 
Plummer external potential. 




GMpi 


(3) 


where Mpi and Opi are the mass and scale length of the Plummer 
potential, respectively. They revealed that if the external potential 
well is deep enough compared to the potential of embedded stel¬ 
lar system, the external potential well behaves like a heat bath. The 
velocity dispersion of the embedded stellar system, therefore, be¬ 
comes isothermal and there exists a quasi-equilibrium solution of a 
potential-density pair for the nearly iso thermal stellar system (see 
Model 2 of Fig. 10 in lYoon et al J2OI ill . 

The bulge is one of galactic components, extending over few 
kpc scales. NC and bulge are independent co mponents of galaxie s 
with different surface brightness profiles (e.g.. fBalcells et alj2003h . 
In the Milky way, the effective radius (i.e., half-light radius) of the 
bulge is about 0.1 kpc, and the mass is estimated to be roughly 
lO^^Mffi. Acco r ding t o a dynamical model of Galactic bulge sug¬ 
gested by I Kent! ( Il992h . the kinematics of the bulge is affected by 
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the MBH, Sgr A*, at inner parsec scale, and the velocity dispersion 
is nearly flat from 1 to 10 parsec and increases gradually at the large 
radii. We assume that the the bulge is a sphere in this study for the 
simplicity. This may not cause a serious effects on the dynamics of 
the nuclear cluster we are considering since the role of the bulge in 
our model is to confine the nuclear cluster within a few parsec. 


3 COMPUTATIONAL METHODS 

In this study, we used th e GPU accelerated version of NBODY6 
(iNitadori & Aarsethll2012ll . This code includes many efficient and 
accurate algorithms such as the fourth-order Hermite integrator, 
the individual and block time steps, the Ahmad-Cohen neighbor 
scheme, the Kustaahe imo-Stiefel (KS) and chain regularization 
scheme iAarsethlll999h . Recently, by using numerous stream pro¬ 
cessors of GPU devices, calculations of gravitational interactions 
among stars have been significantly accelerated through massive 
parallelism. Force calculation in NBODY6 code is divided into regu¬ 
lar force from distant particles whose gravitational potential change 
slowly and irregular force from neighboring particles that give 
strong gravitational fluctuations. In the GPU accelerated version, 
the regular forces are calculated by GPU machines with single pre¬ 
cision while irregular force calculation that needs better accuracy 
is still done by CPUs. 

As we mentioned above, the external potential is composed of 
two parts, 

4>ext = 0MBH -F 4>pl (4) 

where 0 mbh is the Keplerian potential due to the central MBH. In 
our implementation, these external forces are added to the sum of 
regular forces. 

All the calculations in NBODY6 code use dimensionless time, 
length and mass units. From given unit of length f in parsec and 
mean stellar mass M in Mp,, the p hysical unit of velocity and time 
can be expressed as ( lAarsethll2010l) 

_2 /NM 

velocity : 6.557 x 10 ( —^ 



/ ^3 y/2 

time: 14.94 Myr (5) 

where N is total number of stars. 

Table[T]shows model parameters of our simulations. Although 
galactic nuclei contain ~ 10® stars in a cubic parsec, it is hard to 
treat such a large number of particles for NBODY6 code even with 
GPU machine. We, therefore, used several different number of par¬ 
ticles ranging from 10,000 to 100,000 in order to build a scaling 
relation. We also made several simulations with different random 
seeds and took average in order to reduce the statistical noises. 
The masses of MBH are chosen to be 10 and 20 % of the entire 
mass of cluster stars M^i (but excluding the mass of the Plummer 
po tential). Because these mass es are relatively larger than thos e 
of ISigurdsson et al.l h995h and iHollev-Bockelmann et all i2002h . 
we also used longer black hole growth time scale fiviBH = 50. 
The softening parameter of the MBH is fixed to 10“"^, which is 
much smaller than the radius of influence of the MBH (c.f., rinf = 
Mivibh/ctJ ~ 0.1 for Mmbh = 0.1). We assumed that all stars 
have same masses (i.e., m = 1/A^) of stellar mass BHs with the 
mass of IQMq . Thus, the total mass of the cluster in physical unit 



Figure 1. Time evolution of Lagrangian radii for the star cluster with a 
growing central MBH and an external Plummer potential well. By the 
growth of MBH, the Lagrangian radii decrease with time until T = 
Imbh = 50. After full growth of the MBH, the cluster expands due to 
the heating from the MBH and stars in the cusp. 


becomes Mtot = IONMq . It is well known that there is a correla¬ 
tion b etween the mass of MBH s and the kinematics of bulges. Re- 
centlv. lMcConnell & Mai (1201 3h updated the Mmbh — Mbuige rela¬ 
tion and found that the mass of MBH is roughly 0.2% of the bulge 
mass for 10^ < Mmbh/Mq < 10^® (also see I Graham & Scod 
I2OI5I) . For the external potential of a bulge, we fixed the mass and 
scale length of the Plummer potential to 100 and 5, respectively, 
which corresponds to Mmbh /Mbuige = 0.001 and 0.002 for our 
models. Under this potential well, the embedded stellar system is 
expec ted to become nea rly isothermal in a few half-mass relaxation 
time ^Yoon et alj|201 ih . We also considered a model without the 
external Plummer potential well (Model 0) for the purpose of com¬ 
parison. 


4 DYNAMICAL EVOLUTION OE STAR CLUSTERS 
4.1 Cluster expansion 

In order to investigate the effect of the MBH and the surrounding 
bulge, the dynamical evolution of the star cluster is presented in 
this section. Here, we focus on the Model 4 with N = 100, 000 
and Mmbh/Mci = 0.2. 

Fig.Hlshows the time evolution of Lagrangian radii including 
1, 5, 10, 20, 50 and 75 % of total cluster mass. Before T = Imbh = 
50, all Lagrangian radii decrease with growth of the central MBH 
to adjust with the strong potential of the MBH. After T = Imbh, 
Lagrangi an radii increase gradually as reported in previous stud¬ 
ies (e .g., iBaumgardt et al.ll2004al : iFreitag et alJl200d : iFiestas etaP 
I2OI2I) . The basic mechanism of the expansion is similar to that of 
post core collapse expansion of star clusters. Wit h the MBH, ki - 
netic energy can be generated by stars in the cusp dShapirdllPTTh . 
The MBH and the innermost star can behave like hard binaries in 
the core of star clusters. Encountering other single stars in the cusp, 
they are to be bounded stronger and convert their internal energy to 
the kinetic energies of stars in the cusp. The kinetic energies are 
transferred to the whole cluster via relaxation. 

While the kinetic, potential and total energies of an isolated 
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Table 1. Initial parameters for all models. 


Model 

/Vcl 

Nrun 

A/mbh/A/ci 

A/mbh/™* 

Imbh 

£mbh 

Mpi 

flpl 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(V) 

(8) 

0 

100,000 

1 

0.2 

10,000 

50 

10-4 

- 

- 

1 

10,000 

10 

0.2 

2,000 





2 

20,000 

5 

0.2 

4,000 





3 

50,000 

2 

0.2 

10,000 





4 

100,000 

1 

0.2 

20,000 

50 

10-4 

100 

5 

5 

20,000 

5 

0.1 

2,000 





6 

50,000 

2 

0.1 

5,000 





7 

100,000 

1 

0.1 

10,000 






Notes. (1): Number of stars in the nuclear star cluster. (2): Number of simulations with different initial random seeds. (3): MBH mass compared to the total 
mass of cluster. (4) Mass ratio of MBH to the stellar mass. (5): MBH growth time scale. (6): Softening parameter of MBH potential. (7): Mass of external 
Plummer potential. (8): Scale length of external Plummer potential. 



Figure 2. Time evolution of energies for stellar particle only. To be virial- 
ized against the external potential, the kinetic energy becomes much larger 
than the isolated systems. The kinetic energy decreases with time because 
of the cluster expansion while the total energy is nearly conserved. 


self gravitating system are 1/4, -1/2 and -1/4 in NBODY units, 
respectively, the system embedded in a potential well has larger 
kinetic energy tha n the isolated system through the virialization 
dYoon et al.l 201 ill . Also, when the kinetic energy is even larger 
than the magnitude of the potential energy (i.e., the total energy 
becomes positive), the system will not reach the core collapse. Fig. 
1^ shows the time evolution of the energies. Our models are de¬ 
signed to have positive total energy. Initially, the kinetic energy, 
potential energy of self-gravity and total energy are 0.7, -0.45 and 
0.25, respectively. However, as the black hole grows, the poten¬ 
tial energy decreases because the cluster becomes more centrally 
concentrated as shown in Fig. [T] while the kinetic energy increases 
in response to the external potential of the bulge. AtT — /mbh, 
the energies become -0.5 (potential), 0.9 (kinetic) and 0.4 (total), 
respectively. The potential energy increases after T = /mbh due 
to the cluster expansion although we designed a quasi-equilibrium 
model by using the external Plummer potential. Once the MBH has 


been grown, the numerical error is generated by the interaction be¬ 
tween the innermost stars and the MBH at r —>■ 0. For Model 4, the 
root-mean-square (RMS) error of total energy per step (AE/E) rms 
including the potential energies from external potential and MBH 
is ~ 10“®, overall. Such a large error could originate in escaping 
stars which are kicked out from inner region due to the less accu¬ 
rate force calculation during the regularization processes perturbed 
by strong gravitational potential of the MBH. However, they have 
very large escaping velocity (~ lOOa,) and rarely interact with 
other stars during escaping. When these escapers are excluded, the 
RMS error drops to ~ 10~®. Although these escapers rarely inter¬ 
act with other stars, the indirect heating could cause small structural 
changes at r < 10“^ where escapes mostly happen, as discussed 
in S5.2. 


4.2 Radial profiles 


In order to describe the dynamics at the vicinity of MBH, it is very 
important to determine the radius of influence rinf, within which 
stars are directly affected by the gravitational potential of MBH. 
We estimated the radius of influence by following the definition of 
iBaumgardt et alj ( l2004all as 


GMmbh 


( 6 ) 


where ‘2’ is an empirical factor suggested by [Baumgardt et al.l 
(l2004ah . and cr* is here the density-weighted one dimensional ve¬ 
locity dispersion of stars 


(■’’halt 4jYa^{r)p{r)r'^dr 

^ ' _ (7) 

“ 47rp(r)r2dr ^ ^ 

where rhaif is half-mass radius. Since rinf and cr* are interrelated 
each other, the estimation of them is made iteratively within the 
relative error 10”“^. For instance, nnf = 0.16, cr* = 0.79 and the 
mass of stars within rinf, M{< rinf) « 0.3 for Model 4, respec¬ 
tively. 

In Fig. we show the density profile of the star cluster of 
the Model 4 for / = 0 and IOOOTnbody- We see in Figs. [T] and 
l^that the c luster is not in equilibriu m but expanding. However, as 
reported in iBaumgardt et alT(l2004ah . the equilibrium profile is ex¬ 
pected to be established from inner to outer regions after a few lo¬ 
cal relaxation times (c.f., r*! ~ 100 and Trh.o ~ 1600Tnbody for 
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Figure 3. Density profiles at T = 0 and T = IOOOTnbody- Grey line 
shows the initial condition. Because of the existence of MBH, the stellar 
cusp, so-called Bahcall-Wolf cusp, is formed inside the radius of influence. 
Our model has less steep density profile when approaching rinf than the 
theoretical profile of due to the contamination from outer paifs. 


Figure 5. Separation between the MBH and the centre of mass of stars as 
a function of time for simulations with different N. Contrasts show differ¬ 
ent numbers of stars. With less number of stars, the MBH wanders more. 
Dashed lines represent the wondering radius based on the our scaling rela¬ 
tion obtained from W=100,000 model. 



Figure 4. Velocity dispersion profiles at two epochs: T = 0 and 1000. 
The existence of MBH makes stars within radius of influence follow the 
Keplerian profile. For the outer parts, the velocity dispersion is flatter than 
the initial condition due to the external potential well. 


Model 4). The radius of influence Tint is marked as the downward 
arrow. The slope of density cusp at r ^ rinf is similar to -1.75 
for equal-mass systems but becomes shallower at r —>■ rinf. This 
feature has already been observed many previous studies with dif - 
ferent approaches such as dire ct A^-bodv llBaumgardt et al.ll2004^ 
and Monte-Carlo simulations dpreitag et al.ll2006ll . There Ts an up¬ 
turn of the density near r = 5. This radius is the scale length of 
the external Plummer potential, so stars have piled up at this radius 
against the external force from the potential. 

The radial profile of velocity dispersion is presented in Fig.|4] 
The black line shows the profile at T = 1000, and the grey line is 
the initial one. For the outer part in Fig. [d] the velocity dispersion 
is not completely flat, but flatter than the initial profile. There is a 


bump of velocity dispersion at r ~ 5, the scale length of the exter¬ 
nal Plummer potential. This may have been caused by the heating 
effect by the external potential which is more effective where the 
depth of the external potential acts like a reflecting wall. This will 
be more evident in the next subsection on the velocity anisotropy. 
The heating through the reflection would appear mostly in the ra¬ 
dial direction. When the MBH exists, stars inside rinf are strongly 
affected by the MBH, and their velocity dispersion follows the Ke¬ 
plerian profile (i.e., cT{r) ~ However, fhe slope is a lif- 

tle shallower than the expectation as similar to the density profile 
in Fig.|^ Another important characteristic of MBH surrounded by 
stars is the wandering. An MBH in a stellar system can move ran¬ 
domly, like a Brownian motion, by the interaction with stars which 
are bounded and un-bounded to the potential of the MBH. In ad¬ 
dition, the MBH and the innermost star can play a role as binaries 
in the core-collapsed star cluster. They kick a star interacting with 
them with high kinetic ener g y, and this causes the recoil motion of 
the MBH. iLin & Tremaine! (Il980h investigated the motion of the 
MBH in the stellar system such as a globular cluster and concluded 
that the interaction with unboun d stars is the mos t impo rtant effect 
causing the motion of the MBH. fBahcall & WolfI l ll976h estimated 
the uncertainty of the position of the MBH by the wandering of the 
MBH as 

fwand ~ 0.92/‘c yj m.,,/AT mbh (8) 

where is the core radius. It is difficult to define the wandering 
rad ius in our simulatio ns because there is not a well-defined core, 
but iFiestas et alJ (l2012h found the empirical relation for the MBH 
wandering radius as rwand = 0.5(m*/MMBH)'^’* for King mod¬ 
els. 

Fig |5] shows the wandering radius estimated by the radial dis¬ 
tance of the MBH fixed at the origin to the centre of mass of the 
stellar particles. The solid lines are the wandering radius in the 
simulations for the simulations with different number of stars rang¬ 
ing from 20,000 (top) to 100,000 (bottom) with the MBH mass of 
0.2. As the number of stars becomes larger, the wandering radius 
becomes smaller as shown in equation The dashed lines are 
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Figure 6. Time evolution of the slope of stellar cusp within the radius of 
influence (see text). Different lines show results with different numbers of 
stars. There is a tendency of steeper slope for the larger N models. The 
model with A'^=100,000 reaches the Bahcall-Wolf slope at T 2000. 



0.01 0.10 1.00 10.00 


the values obtained with the relation from iFiestas et alj ( |2012|) . but 
with different coefficient of 0.7 rather than 0.5. Thus the wondering 
radius follows the scaling relation with N very well. Different co¬ 
efficient is likely to be from the different initial prohles. The time 
evolution of the slope of the central stellar cusp 7 is shown in Fig. 
| 6 ] Due to the fact th at the density cusp is smeared out by the ran¬ 
dom walk of MBH jFiestas et al.ll20la) and the density profile is 
contaminated by outer region at r —>■ rinf as shown in Fig. [ 3 ] the 
slope is estimated with the stars at rwand < r < O.Srinf (However, 
we used larger factor rather than 0.3 for small N runs to secure 
enough number of stars in that range, and this may cause shallower 
cusp for small N runs). Different contrasts mean different number 
of stars. For larger N, the slope increases with time slowly toward 
the Bahcall-Wolf value of -1.75 while those with smaller N do not 
get close to the Bahcall-Wolf cusp. For Model 4, the cusp once ap¬ 
proaches to the Bahcall-Wol f cusp at T = 2000rNBODY. This is 
consistent with the result of [Freitas et alj ( l200d) that it takes or¬ 
der lOrri (c.f., Tri ~ 100) for the Bahcall-Wolf cusp to be fully- 
developed. However, it is still true that most of our models have 
shallower cusps overall. The shallower stellar cusp will affect the 
merger rate of stellar mass BHs near the MBH as discussed in the 
next sections. 


4.3 Velocity anisotropy 


It is well known that the radial anisotropy in the velocity dispersion 
increases at the outer part of isolated stellar systems as a result 
of two-body relaxation dOiersz & Heggielll9 ^:ISpitz eJl98'i|). The 

ain^ ' 


anisotropy parameter can 


be defined by 1 Binnev & Tremain 320081) 


P 



(9) 


where crt and at are the tangential and radial velocity dispersion, 
respectively. This anisotropy parameter varies from —00 (purely 
circular orbits) to 4-1 (purely radial orbits). Also, in a tidal field, 
the radial anisotropy decreases during post core-colla pse expansion 
due to the loss of radial orbits dTakahashi et aljfl997h . However, in 
the case of stellar systems with a growing central massive object. 


Figure 7. Profiles of Radial and tangential velocity dispersion (upper) and 
the anisotropy parameter (lower) for Model 0 without the external Plummer 
potential. The velocity distribution is ta ngentially biased, an d the profile is 
similar to that of the isochrone model in lOuinlan et al.l d 19951) except for the 
position of maximum anisotropy. 


the tangential anisotropy i s developed (i.e ., P < O.lYoundl 19801: 
Goodman & Binned 1 19841: [O uinlan e t al.l 1 199.51: ISigurdsson et al.l 


19951 ; iHollev-Bockelmann et al] 2002l) T Quinlan et al. dl995t) have 

revealed that the aspect of anisotropy is affected by the initial 
models. For models with a core such as isothermal sphere and 
isochrone model, the velocity distribution become s isotropic in 
the limit r —)■ 0 dYoundll98(ii : lOuinlan et alJll99^ . On the other 
han d, for ‘two-pow er’ models for galaxies like Dehnen’s mod¬ 
els dPehnenl I l993h . the velocity distribution is still tangentially 
biased at r —» 0 dOuinlan et al.l 1 19951 : ISigurdsson et al.l 11993 : 
iHollev-Bockelmann et alj2002l) . The tangential anisotropy reaches 
the peak at r ~ rinf, and the veloci ty distribution becom es isotropic 
at r 3> finf again (see Figs. 2-5 of lOuinlan et alj 19951) . 

The radial and tangential velocity dispersion and the velocity 
anisotropy parameter after the growth of MBH are shown in Fig. [7] 
for the Model 0 without the external Plummer potential and Fig [8] 
for the Model 4 with the potential, respectively. For the Model 0, 
the velocity dispersion decrease rapidly with radius unlike those of 
the model with the potential well. Although we use the Plummer 
model as the initial model, the anisotropy parameter is likely to be 
simila r to that of the isochrone model as shown in lOuinlan et al.l 
)ll995h . However, the maximum anisotropy is located at the larger 
radius than rinf. On the other hand, for the Model 4, the radial ve¬ 
locity dispersion is enhanced due to the radial acceleration from 
the external potential as mentioned in the previous section. Obvi¬ 
ously this is an artifact of the fixed external potential. If the external 
potential is a live one, the velocity dispersion and the anisotropy 
would be more smooth. However, as we will see in the following 
section, the details of the velocity profiles and anisotropy would not 
affect on the estimation of binaries as the binary formation rate at 
the outer part becomes very low. 
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Figure 8. Profiles of Radial and tangential velocity dispersion (upper) and 
the anisotropy parameter (lower) for Model 4 with the external Plummer 
potential. For the inner region, the velocity distribution is similar to that of 
isolated models. However, the tangential anisotropy is driven by the external 
potential well at the outer region. 


5 BLACK HOLE BINAIHES 

5.1 Close encounters and GR capture 


In order to estimate the merger and detection rates of BH-BH bi¬ 
nary coalescences, we need to know binary formation rates as well 
as the orbital parameter distribution just after the capture. We only 
consider the direct capture by gravitational radiation during close 
encounters. 

The energy loss and change of orbits by GW f or binary sys¬ 
tems were first studied by IPeters & MathewsI ( 19 ^ ) b ased o n the 
post-Newtonian ( PN) approximation. Later , iHanse 3 (Hzl) ex- 
tended the study of IPeters & Mathe^ ( Il963h to the hyperbolic en¬ 
counters. With given masses mi, m 2 , a semi-major axis a (defined 
as a = Gm\m 2 l‘lEQ where Eq is the initial orbital energy) and 
an eccentricity e, the energy and orbital angular momentum losses 
by GR are given by 


AE = 


X 


2 mfm2(mi-b m2)^^^ 

~T5~i? a7/2(e2 _ 1)7/2 

(tt — 0o)(96 + 292e^ -I- 37e^) + iesin0o(6O2 + 457e^) , 

( 10 ) 


AL, 


o 2 2 

8 G mim2 

5 c® o2(e2 - 1)2 


(tt— 0o) (8-l-7e^)-l-e sin 9o{13+e^) 

( 11 ) 


where G, c and 9o are gravitational constant, speed of light and the 
incidence angle at infinity defined as do = cos“^(l/e), respec¬ 
tively. Two encountering but unbound stars, therefore, become a 
binary if the energy loss by GR is larger than the orbital energy Eq. 


From equations Gil and GB, one can obtain the semi-major axis 
and the eccentricity of the captured binary as 

, Gmim2 


2{Eo + AE) 


e' = Wl- 


(mi -b m 2 ){Lz,o + AL^) 


Gm\r 


( 12 ) 


(13) 


where the subscript 0 indicates the initial value. 

Many previous studies have discussed t he GR captures of 
compact stars in dense stellar s ystems (e.g., lOuinlan & Shapir^ 
Il98ill989l: l0’ Leary et ^ |2009|) . The stating p oint is the cross 
section for GR capture. Quinlan & Shapircl ( ll987l) deduced the cap¬ 
ture cross section under the parabolic approximation. This approx¬ 
imation is valid because trajectories of the stars near the pericentre, 
where bulk of the GWs is radiated, are almost identical to parabolic 
with the same pericentre distance. The equation llOl l is rewritten 
with parabolic approximation as 


AE = - 


SStt G’’^^ mim|(mi -b m2)^^^ 
12^2 c® 


•J12 


(14) 


where rp is the pericentre distance. Again, the GR capture will hap¬ 
pen when the energy loss is lager than the orbit al energy. By the re¬ 
quirem ent of I AD I > I2{m\ -bm^l. lOuinlan & Shapir3 

J 19891) obtained the maximum pericentre distance for GR capture: 


857rx/2 G’^/^ mim2(mi -b 


12 


2/7 


(15) 


where v^a is the relative velocity at infinity. Therefore, the capture 
cross section is given as 

2G(mi -b m 2 ) 


“^cap — ^^p,max 


' p,maxL'oo 


- 17 


G^niimori ®^’^ 


(16) 


where rj is the symmetric mass ratio defined as t) = mxm^l{jn\ -b 
mo)^■ We assumed that gravitational focusing is dominant com¬ 
pared to the geometrical cross section for the last equality in the 
above equation. 

Here, we are going to introduce a statistical interpretation of 
GR captures to understand the situations and predict the BH-BH bi¬ 
nary coalescences in realistic regime. We can assume that the mo¬ 
tions of stars follow the one-dimensional normal distribution with 
a given velocity dispersion a. From the equation l ll6t , the distribu¬ 
tion of the pericentre distance of encountering stars also becomes 
uniform (i.e., dS = d{nb^) ~ 27rG(mi -b m 2 )/w^ • dvp) if the 
gravitational focusing dominates. Therefore, for unbound close en¬ 
counters that lead to the formation of binaries by GR, rp/rp.max 
is distributed uniformly in the range [0, 1]. Under the parabolic ap¬ 
proximation, the semi-major axis and the eccentricity can be rewrit¬ 
ten with a and rp/rp.max from the equations 1141 and IIS! as 


GM 


1.58AU 


M 


2 OM 0 


75km/s 


(17) 


e = 1 -b 


f 170777)^ 


2/7 


X 3 ) 

1 -b 1.57 X 10“®?)^^'^ 


10/7 


10/7 


75km/s 


( 18 ) 
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Figure 9. Distribution of semi-major axis (left panel), eccentricity (middle 
panel) and merging time (right panel) for WMq BH-BH binaries in the 
Milky-Way-like galaxies (i.e., a = 15 km/s). 


where M is the sum of masses, and we set Uinf = \/2a. By assum¬ 
ing rp/rp,max = 1/2 typical pericentre distance for GR capture of 
encountering two IOMq BHs is 2514km, corresponding to about 
840 -Rsch (Schwarzschild radius), in Milky-Way-like galaxies (i.e., 
(T = 75km/s). Because the pericentre distance is nearly the same 
before and after capture, the semi-major axis and eccentricity of the 
binary formed by GR capture are given by 


- 7/2 


- 1 


0.153AU 


M 


20Mq ) \ 75km/s 


1 - 1.62 X 


75km/s 


10/7 


(19) 


( 20 ) 


For the case with velocity dispersion of cr ~ 75 km/s, the eccen¬ 
tricity of a typical binary formed by GR capture is 1 — e' ~ 10“"^. 
From the distribution of semi-major axis and eccentricity, we can 
determine the distributio n of merging time for such binaries. The 
merging time is given by jPeterJ 19641) 

73 2 37 4 I 


Tmer = 


5 

64 G®mim2(mi + m 2 ) 




( 21 ) 


where ao and eo are the initial semi-major axis and eccentricity, 
respectively. Since the orbits of binaries in our consideration are 
nearly parabolic, the merging time has a strong dependence on the 
velocity dispersion, i.e., Tmer ~ ao(l — ~ Note 

that this merging time is not corrected for lower order PN terms 
(see §7 and Appendix A for other PN corrections). In Fig. we 
have shown the distribution of semi-major axis, eccentricity and 
merging time of BFl-BH binaries from equations ( I19t , OOb and 
(ED. As mentioned, we assume that the velocity of star follows 
one-dimensional normal distribution with a= 75 km/s (i.e., the ve¬ 
locity dispersion for the Milky Way) and the pericentre distance 
t'p/i'p.max follows Uniform distribution. The peak position of each 
distribution is equivalent to the typical value in equations l ll9b and 
(Ho}. For merging time, the mode is ~ 10® second. Although the 


hyperbolic encounters 



10“® 


10“" 10“® 10“ 10® 

eccentricity (e—1) 


Figure 10. Distribution of semi-major axis and eccentricity of GR captured 
binaries for the Model 4. Different colors mean the velocity dispersion at 
the radius of influence in physical units as indicated by the color bar. These 
encountering stars become binaries by GR capture. 


distribution is quite wide, almost all binaries will merge within a 
few thousand years. 

5.2 Merger rates 

We collect the parameters of all close encounter events in our sim¬ 
ulations in order to investigate the GR capture and the compact 
binary coalescences. In the NBODY6 code, the regularization al¬ 
gorithm helps the calculation of very close orbits with high pre¬ 
cision. If the separation or the time step of stars becomes smaller 
than certain criteria, the stars are taken away from the main loop, 
and their motions are calculated with time smoothing. We turn on 
the KS regularization, the two-body regularization scheme, and ex¬ 
tract the semi-major axes and the eccentricities of pairs experienc¬ 
ing close encounters at the pericentre passage to avoid the effect 
of perturbation by nearby stars. We have shown an example of 
(a, e) distribution of the close hyperbolic encounters for the Model 
4 after T > /mbh in Fig[T0] Each filled dot represents each en¬ 
counter and the dashed line shows the limit of KS regularization 
(i.e., Tperi ~ 10”“^), SO pairs lying above this line are not consid¬ 
ered in our investigation. In order to determine whether a certain 
encounter results in a binary, we therefore need to convert our di¬ 
mensionless results to physical quantities according to the equation 
(EJ. In Fig. [To] there are some orbits with rainbow colors. These 
colored orbits will become binaries when the overall velocity dis¬ 
persion of stellar system is larger than that velocity (e.g., the orbit 
colored red at the left-end will become a binary when the veloc¬ 
ity dispersion is ~400 km/s.). By counting the number of capture 
events, we estimate the binary formation rates in our simulations. 

For a single BH passing through stars with a speed v, the time 
scale for GR capture is 

/cap = (Ecapnu)”^ (22) 

where n is the number density of background stars. The binary for¬ 
mation rate between stars with different mass mi and m 2 in the 
shell with the range [r, r + dr] can be expressed as 

= 47rr^ni(r)n2(r)(Ecap«) 


(23) 
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log r 


Figure 11. Radial distribution of GR binary capture rates. The solid curve is 
the semi-analytic model from the equation (25). This semi-analytic model 
agrees well with events in the simulation except for the inner region. From 
the theoretical models, the slope of distribution inside the radius of influence 
is expected to be 2/7. Although the actual slope differs from the theoretical 
models due to the overall incompleteness of Bahcall-Wolf cusp as shown in 
Fig.[6] the contribution to the overall merger rate within this radius is small. 
Thus, it does not significantly affect to the estimation of the entire merger 
rates. 


where (Ecapi^) is the velocity averaged value. Thus, assuming 
t^oo = y/2v in equation GD and replacing v by (j(r), the veloc¬ 
ity dispersion, we have 

~ 87--r ni(r)n 2 (r)cr(r) (24) 

For the case of systems composed identical stars, m = mi = m 2 , 
this equation becomes 


(25) 


where the half is to avoid double counting and 77 = 1/4. Fig. 
m shows the radial distribution of cumulative capture events 
dNcap/dlogr during At = 2000Tnbody- Flistogram is from N- 
body simulations of Model 4, and the noisy line is from the equa¬ 
tion IB) (i.e., rfVcap/dlogr = rdrcap/dr ■ At) with the density 
and the velocity dispersion profiles from Figs.j^andj?] respectively. 
In order to get large sample size, we set the unit of velocity to the 
half of the speed of light (0.5c). Although this velocity is unrealis¬ 
tically high, it is possible to guess what happens in realistic situa¬ 
tions because there is no relativistic effect on the simulations. The 
simulation results and our semi-analytic model from given density 
and velocity dispersion profiles show good agreement at the radius 
larger than rinf. However, at r < rinf, there is some discrepancy 
between them: the binary formation rate obtained with simulation 
is smaller than that with the equation l l25t . One possible reason is 
the time variation of density structure within the radius of influ¬ 
ence as shown in Fig. From the theoretical model of star dis¬ 
tribution within the radius of influence (e.g., p(r) ~ 


cr(r) 


r-1/2 


), we obtain the slope of dNcap/d log r as 


dN, 


cap 2/7 

oc r ' . 


. (26) 

dlogr 

However, the slope inside rinf in Fig. [TT]is not the same as 2 be¬ 
cause of the incompleteness of the Bahcall-Wolf cusp (see Fig. 


[ 6 ). Around the wandering radius, even merger events from N- 
body simulations and the semi-analytic model show disagreement. 
This gap is likely due to escapers that cause indirect heating at 
that radius. Incidentally, more than 80 per cent of events occurred 
outside Hnf, and the peak of dNcap/dlogr is located at rhaif- 
The fact that the most of the merger takes place outside the ra¬ 
dius of influence is mainly because the density distribution fol- 
lows Bahcall-Wolf or shallower. Our results are quite different from 
lO’Learv et al.l (20091) who found that the merger rate is more cen¬ 
trally concentrated because the distribution of the black holes is 
much ste eper than Bahcall-W olf due to the mass segregation (see 
Fig. 7 of lO’Learv et al.ll2009ll . Our simulation is one component, 
and thus cannot reproduce the effect of the mass segregation. We 
note here that many recent simulations with different approaches 
showed that the density profile of the massive component fol¬ 
lows the Bahcall-Wolf or steeper while lower mas s component fol¬ 
lows sh allower slope (direct A^-bod y simulations , Baumgardt et al.l 
l2004bh (Monte-Carlo simulations, [Freitag et al.ll2006lf and (time- 


dependent Fokker-Planck equation. Fiestas et al., in preparation). 
If that is the case, the small discrepancy of analytical estimation at 
small radii thus does not affect the estimation of the total capture 
rates. 

In order to obtain the overall merger rate for NC, we need 
to integrate the equation I25t over the volume. We can assume that 
the merger rate is equivalent to the capture rate because the merging 
time of BH-BH binary in our simulations is negligible compared to 
the cluster time scales as discussed before. It is difficult to estimate 
the merger rate because n{r) and (j{r) are not simple function of r. 
Assuming that the velocity dispersion remains a constant over the 
entire cluster, and replacing the integral J rdr^dr by An, we can 
write the integrated merger rate as 


Finer ~ ■ N ■ u ■ a^. 

~ (27) 

where n, cr* and M are mean number density, the velocity disper¬ 
sion of the system and total mass of the cluster, respectively. We 
have used the virial theorem to get the last relationship. We see 
that the event rate is inversely proportional to the total mass of the 
cluster with rather steep dependence on the velocity dispersion of 
the cluster. To convert our results to physical units, it is necessary 
to determine the representative value of velocity dispersion of A- 
body simulations. We already estimated the density-weighted ve¬ 
locity dispersion in equation 0 as simila r to the observational esti - 
mation of velocity dispersion of bulge as (McConnell & Ma|l2013h . 
Note that the effective radius of bulge that is the upper limit of the 
integration for the systemic velocity dispersion in observations is 
much larger than the half-mass radius of NCs. According to re¬ 
cent obse rvation for the cal ibration of velocity dispersion of nearby 
galaxies (Kang et^l2013h . however, the velocity dispersion does 
not change much with the aperture size. Thus, the velocity disper¬ 
sion of NCs is nearly identical to that of bulges. 

Figs. [T^ and [T^ show the merger rates as a function of the 
velocity dispersion in the physical unit. The mass ratio of MBH 
to the cluster is fixed at 0.2 for Fig. [T^fcr* ~ 0.79) and 0.1 for 
Fig.[T3](a* ~ 0 .75), respectively. Different symbols represent the 
different models with different A and dashed lines are from the 
time averaged result of numerical integration of the equation (27). 
Because of the limitations in the number of particles and the inte¬ 
gration time, we only consider the unrealistic range of the velocity 
dispersion (l,000km/s < cr* < 20,000km/s). Nevertheless, because 
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Figure 12. Merger rates as a function of velocity dispersion for models 2- 
4. The MBH mass is 20% of the total mass of the cluster. Filled symbols 
are from the number counts in simulations with different number of stars. 
Dashed lines show the equation GiJ with the proportional constants ob¬ 
tained from the simulations. There are good correlations between merger 
rates and velocity dispersion. 



t/^rl,0 


Figure 14. Time evolution of merger rate from Model 2 as a function of 
time at the radius of influence. The merger rate is scaled for a Milky-Way¬ 
like galaxy. Noisy line and histogram are from the integration of equation 
IID and the number of events counted in the simulations, respectively. Due 
to the cluster expansion, the merger rate decreases with time. The dashed 
horizontal line represents the merger rate from equation {28). 



Figure 13. Merger rates as a function of velocity dispersion for models 5-7. 
The MBH mass is 10% of the total mass of the cluster. 


the results show very good scaling relation with the velocity disper¬ 
sion and the total mass of the cluster as predicted by the equation 
(|27j, it is possible to extrapolate our results to realistic parameters 
for NC. As a result, the merger rate can be expressed by an equation 
with the mass of MBH and the velocity dispersion: 

-1 


^2.06 X 10”‘‘Myr"^ 


A/mbh 


Y. 


3.5 X 106 M© J I 75km /s I 


\ 31/7 


Thus, the merger rate is about 2.06x 10“^°yr“^ for Milky-Way¬ 
like galaxies if we assume that the total mass of embedded star 
cluster is 5 times heavier than the MBH. 

Our realization of the NC is not static: the NCs expand slowly 
with time even though they are bounded by the external potential as 
described in §4. The time evolution of the merger rate for a Milky- 
Way-like galaxy is represented in Fig. [14] which is estimated from 
the Model 2 with N = 20, 000 in order to see the long-term evo¬ 
lution. The time i s scaled by the initial relaxation time at the radius 
of influence Tri,o ( ISpitzeJl987h after growth of the MBH 

Ib.dG^mpinf In A ^ 

where fh and pinf are the mean mass (for equal-mass, fh = M/N) 
and the mean density inside the radius of influence, respectively. 
The noisy line is the merger rate from the numerical integration of 
the equation ( 125b . For comparison, we plot the histogram showing 
the number of events counted in the simulation. They are selected 
with (T* = 10"^ km/s for sufficient samples and rescaled to ct* = 75 
km/s range. The horizontal dashed line is the time-averaged merger 
rate in equation l l28t . Due to the expansion of the cluster, the merger 
rate decreases with time. These two estimates show good agree¬ 
ment, and therefore it is possible to surmise the merger rates from 
given density and velocity structures of stellar systems. The conver¬ 
sion of these merger rates to the detection rates for GW detectors 
will be presented in next section. 


6 DETECTION RATES OE BH BINARY MERGERS 


for Mmbh = 0.2Mtot, 


(28) 


Finer «8.58 X 10"®Myr"^ 


Mmbh 
3.5 X 106 Mq 


-1 


(T* 1 

75km/s j 


31/7 


(29) 


To determine the detection rate of GWs from BH-BH binary coa¬ 
lescences for GW detectors, it is necessary to calculate how many 
events occur per unit volume in the universe and horizon distance 
of GW detectors. In the previous section, we estimated the merger 
rates in NC as a function of the mass of MBH and the velocity dis¬ 
persion. It is well known that there is a good correlation between 


for Mmbh = O.lMtot. 
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the mass of MBH and the velocity dispersion of surrounding stars 
('e.g. jTremaine et alj2002h 

Mmbh « 1.3 X lO®M0(cr*/2OOkms“^)'‘ (31) 

in a range of the mass of MBH [10 ®Mq, lO^M©]. Later, 
iBarth et ^ bOOSi) confirmed that the relation is also valid for 
MBHs down to 10® Mq. The merger rate for a NC, therefore, has 
a weak dependence on the mass of central MBH T ~ .^mbh 
dO’Learv et alj2009l) . and the merger rates of l l28b and ( 129b become 


Tmer « 3.33 X 10"'‘Myr"^ 
for Mmbh = 0.2Mtot, 


Mmbh 


3.5 X 10®M0 


3/28 


(32) 


1.39 X 10"'‘Myr“^ 


Miv 


3.5 X 10®M0 


3/28 


for Mmbh = O.lMtot. 


(33) 


However, there are several factors that give rise to uncertain¬ 
ties in merger rates. From the equation l l24b . we can infer that the 
merger rate is proportional to the total mass of different mass com¬ 
ponents Ml and M2. As we mentioned before, however, we as¬ 
sumed that all stars are IOM0 BHs. There exist other stellar ob¬ 
jects such as M S stars, white dwarfs (WDs), NSs and BHs in real 
stellar systems. iHopman & Alexanden ( I2OO6I) have studied the ef¬ 
fect of mass segregation of stars around a MBH and concluded that 
the number fraction of different stellar objects evolves from the 
initial state (i.e., A^ms '• (Vwd '• A^ns '• Abh = 1 : 0.1 : 0.01 : 
10“®, for continuously star-forming populations: I Alexandeill200^ 
to Ams : Awd : Ans : Abh = 1 : 0.09 : 0.012 : 0.06 within 0.1 
parsec for Milky-Way-like galaxies. If we set the mass of MS stars 
(O.7M0), WDs (O.6M0) and NSs (I.4M0), the mass fraction of 
BHs A4bh is about 44 per cent of the total mass. When we simply 
assume that the merger rate of BHs in galactic nuclei can be ex¬ 
pressed by the equations l l32b and l l33b with multiplication of AIbh^ 
the merger rate in equation l l32b is reduced to 6.4 5x 10~^^yr~^, 
which is about 3 times smaller than the estimation of lO’Learv et al.l 
(l2009h for BHs with similar mass ranges. Of course, it is more 
complicated to correct for the mass function rather than our con¬ 
sideration because the mass fraction of BHs varies with the radius. 
Furthermore, the mass fraction in our consideration is adequate 
for innermost region although the capture events happen most fre¬ 
quently around the half-mass radii as shown in Fig. [TT] The mass 
fraction of BHs around the half-mass r adius might be smaller than 
that from iHopman & Alexandej ( l2006h . and thus, our results could 
be an overestimation. 

As discussed earlier, our model is one component, and the ef¬ 
fect of the density profiles around the central black hole may not be 
correct. If the density profile of the massive component is steeper 
signific antly for massive com ponents, t he merger rate could b e en¬ 
hanced dO’Learv et alj200^ . Moreover lO’Learv et al.l ( l2009ll have 
shown that using different mass spectra of stellar mass BHs can 
affect not only the merger rates by a factor of ~2 but also the de¬ 
tection rates by a factor of ~10. In that sense our result could be 
an underestimation. However, as we mentioned. Fiestas et al. (in 
preparation) found that the density profile of BHs does not signifi¬ 
cantly deviate from that of dominant components (low mass stars) 
from time-dependent Fokker-Planck simulations. 

The dynamical evolution of NCs also affects the merger rates. 
The merger rate varies at most by a factor of ~2 from T = 0 to 


T = 200rri,o as shown in Fig. [14] iMerritt et~^ (|2007|) have esti¬ 
mated the re laxation times Tri .o for ACS Virgo samples of galaxies 
observed bv ICote et al.l (l2004l) . and found the relation between the 
relaxation time and the central velocity dispersion. According to the 
relation, rri,o is less than a Hubble time with smaller velocity dis¬ 
persion than 100 km/s, corresponding to Mmbh ~ 1.6 x 1O^M0. 
Therefore, the merger rates for NCs with smaller MBHs can be af¬ 
fected by the dynamical evolution. In addition, the relaxation time 
of galactic nuclei implies that most of NCs with larger MBHs do 
not contribute to the merger rates as much as those with smaller 
MBHs because the number fraction of BHs in relaxed nuclei is 
several tens of times lager than that of initial co nditions due to 
the mass segregation dHopman & Alexandej200^ . and the merger 
rate weakly depends o n the mass of central MBH Jo’Learv et al.l 
l2009h . lO’Learv et alj (l2009l) also noted that the variance of the 
number density of galactic nuclei can affect the merger rate. They 
ha ve estimated the var iance of the number density from the results 
of IMerritt et al.l ( l2007l) and found that the merger rate is e nlarged 
as mu ch with the rescale factor ^ ~ 10 — 100. However, see lTsand 
ll2013l) . for possible reduction of the rescale factor by about a factor 
of 5. 

In order to calculate the merger rate per unit cosmological 
volume, we convolve the merger rate per NC with the number 
de nsity of MBHs in the universe (for more de tails, see §3.3.5 
of lO’Learv et aljfioogll . lAller & Richstonel ( l2002h determined the 
number density of MBHs from the luminosity function of galaxies 
as 

^^mbh ^ ( Mmbh \ Mmbh/m. 

dMMBH M. ) 

with the best fitting parameters of = (3.2 x 

10~“M~^Mpc ~®, 1.3x 10® M0,1 25). N ote that using the MBH 
mass function of lAller & Richstona ( l2002h can lead to an overesti¬ 
mation of the detection rates at least by a factor of 2, especially at 
the low mass range of MBHs dOraham & Drivej2007l) . The merger 
rate per volume, therefore, is obtained by integrating the rate over 
the SMBH mass distribution 

TZmeT = f rmer(MMBH)C-7T7—— rfil^MBH (35) 

J Ml dlVlMBH 

where Mu and Mi are the upper and lower limits for integration, 
and Fm.gai and ^ are the merger rate per galaxy as a function of 
the mass of MBH and the scale factor for number density variance, 
respectively. The upper limit can be fixed to Mmbh ~ lO^M© 
by the time scale requirement as discussed above. However, we 
still do not have the exact lower limit of the MBH mass, which i s 
currently about 10® M0 from the observation of iBarth et al.ld2005h . 
By taking the lower limit of the MBH mass to 10"^M0, the equation 
l l35t gives us the merger rate density 

Turner ~ d.lFmer.MwCaoMpC ® (36) 

where Fmer.MW is the merger rate for a Milky-Way-like galaxy, 
and ^30 is the rescale factor for the variance of the n umber density 
of sta rs normalized by 30 (i.e., 1/3 ^ fso ^ 10/3: IO’Learv et aH 
I2OO9I1 . Choosing different lower limit of the MBH mass from 
10® M0 to 10® M0 will give us an uncertainty by a factor of ~3 
in the merger rate density. 

Now, we can estimate the detection rate of BH-BH binary 
coalescences by next generation GW detectors. By assuming that 
the merger events occur uniformly in the universe, the detec¬ 
tion rate only depends on the size of cosmo logical volume which 
we can cover and can be expressed by dO’Learv et al.l 1 20061 : 
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(Belczvnski et alJUoOTl : lOowning et alj20ll] ; ISae et alj2014h 

(37) 

J 1 + z dz 

where z is the cosmological redshift, and the factor of (1 + z)~^ 
represents the cosmological time dilation. For existing GW detec¬ 
tors, the effect of redshift can he negligible because their cover¬ 
age is not too far (i.e., the horizon distance Dh are 33 Mpc for 
NS-NS binaries and 161 Mpc for BH-BFl binaries correspond¬ 
ing to z ~ 0.01 and 0.04 in standard ACDM cosmology, respec¬ 
tively; lAbadie et aLlbOld) . However, for next generation GW de¬ 
tectors, the effect of redshift becomes important, especially for BH- 
BH binaries. The maximum horizon distance Dh can be obtained 
from signa l-to-noise ratio (SNR ) of GW signals (for more details, 
see §4.2 of lO’Learv et al.ll2009l) . Because the redshift affects both 
the mass of source and the frequency, SNR sh ould be estimated 
from the waveforms carefully. Only few st udies dBaker et al.l2007l : 
lO’Learv et al.ll200^ ; iReisswig et al.ll2009li have estimated Dh for 
BH-BH binaries for given SNR. 

Table shows the detection rates expected for advanced 
LIGO. Because we mainly considered the equal-mass models for 
only BHs, the detection rates are corrected by the factor of A4 bh- 
We estimated detection rates by using Dj, from different st udies 
teaker et aD l2007l: [Reisswig et al.|[20OTl; lAbadie et al.l|2010|) . All 
DhS are corrected f or the orientation of sources. Note that Dh from 
lAbadie et alJ ( l201C ) does not include the cosmological effect of the 
redshift. Dh from Reisswig et alJ (l2009h is for spinning BHs but 
independent of the spin of BHs for IOMq BHs. We list three de¬ 
tection rates 7?.det,i, T^-det.re and 7?.det,h in Table for different 
models and different estimates of detection horizon (Dh). We first 
calculated the ‘reasonable’ expected detection rate 7?.det,re using 
the merger rate given in equation l l36b . Given several uncertainties 
including the variance of the number density, dynamical evolution 
and the lower limits of the MBH, we assume that the total range 
of uncertainty is a factor of 60. Thus ‘low’ expected rate 7?.det,i 
is simply obtained by dividing 7?.det,re by x/M while the ‘high’ 
expected rate 7?.det,h is calculated by multiplying the same factor. 
The whole range of the expected detection rate lies between 0.02 
~ 14 yr“^ depending on the maximum horizon dista nce and uncer¬ 
taintie s. These estimates are able to cover those of lO’Learv et al.l 
(I2OO9I) (5-20 yr“^) for BHs around IOM0 although the lowest 
value is quite smaller than that. 

Our estimations have some limitations; (1) We ignore the ini¬ 
tial mass function. This mass function may affect not only the evo¬ 
lution of systems by the relaxation between mass components but 
also the merger rates for BHs with different m asses. (2) We need 
to co nsider various range of Mmsh/Mci (e.g.. lGraham & Snitleil 
I2OO9I suggested that there is a rough relation between the mass of 
MBH and NC.) (3) The exact calculation for SNR is necessary in 
order to obtain more reasonable detection rates. We considered the 
mass segregation and its effect only approximately. These limita¬ 
tions would be considered in future works. 


7 BLACK HOLE BINARY COALESCENCE AND 
WAVEFORM 

Solving Einstein field equation exactly is very difficult and has only 
been done numerically. Fortunately, during inspiral phase of com¬ 
pact binary coalescences, the Einstein equation can be simplified 
with the PN expansion. When a compact binary is formed, the or¬ 
bit decays with time due to the GR. The orbit-averaged change of 



0 5 10 15 20 25 30 

Time (min) 


Figure 15. Time evolution of semi-major axis and eccentricity of a coalesc¬ 
ing binary. Solid line is from the integration of Peters formula. Open and 
filled circles are the results of simulations with 2.5PN correction only and 
full PN correction, respectively. The simulation with 2.5PN correction only 
agrees well with Peters formula while that with full PN correction signifi¬ 
cantly festinates compared to Peters formula. 


the semi-maj or axis and the eccentricity by GR is first derived by 
IPetersId 19641) as 

/ da\ 64 -I-m2) / ,73 2 ,37 4\ 

c5a3(l-e^)r/2 [^+^ 4 ^ 

/ de\ 304 G^mim 2 {mi + m 2 ) f 121 2\ 

\^/ c5a4(l - e2)5/2 + 

in the 2.5PN order (~ 1/c®, the first order GR term). However, 
the orbital evolution is also affected by other PN order terms 
such as IPN (relativistic precession), 1.5PN (spin-orbit coupling), 
2PN (spin-spin coupling, high order relativistic precession) and 
hig her orders. With full c onsideration of PN terms up to 2.5 or¬ 
der, [BMentzen^eT^ ) l2009l) noted that the decay of the binary orbit 
is much faster than that with 2.5PN only. Although the effect of the 
spin is quite important for the motions and waveforms of BH-BH 
binary coalescences, we only consider non-spinning BHs and take 
1, 2 and 2.5PN order terms in this study. 

The equation of motion with PN correction up to 2.5PN 
order can be simply written in the cen tre of mass frame as 
jBlanchet & Iveill2003l:lMora & Willl2doi) 

Gm Gm , , r „ , 

a = an-l-apn =-—r H- —{A-+Bw) (40) 

where A, B are PN coefficients depending on their masses, the rel¬ 
ative position r and relative velocity v (see Appendix A, for more 
details). Many authors have incorporated the PN corr ected force 
in direct W-body simulations with d i fferent PN orders jLeell 19931 ; 
I AarsethlboOVl : l^rentzen et al.ll200j : lBrem et alJl2014h . Similarly, 
we implemented the PN equation of motion to the KS regulariza¬ 
tion process in NBODY6 code. In the KS regularization process, 
a two-body motion is sometimes perturbed by other neighboring 
stars, and these perturbation should be corrected. Thus, we can con¬ 
sider the PN force as a perturbing force in the code by adding the 
PN force and its time derivative. We designed that the binary will 
merge when the separation is smaller than four Schwarzschild radii 
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Table 2. Detection rates of BH-BH binaries for advanced LIGO. 


Models 

Mmbh/Mci 

Fmer,MW 

(Myr-1) 

^BH 

Uh 

(Mpc) 

^det,l 

(yr-h 

''det,re 

(yr-l) 

^det h 
(yr-l) 





986*= 

0.05 

0.39 

3.0 

1-4 

0.2 

3.3 X lO-'^ 

0.44 

~1100^ 

0.07 

0.51 

4.0 





~19009 

0.23 

1.8 

14 





986*= 

0.02 

0.16 

1.3 

5-7 

0.1 

1.4 X 10““^ 

0.44 

~I100l' 

0.03 

0.21 

1.7 





~19009 

0.10 

0.74 

5.7 


“Mean mass fraction of BHs from lHopman & AlexandeJ i2Q06l) . 

6,c,dLow, realistic and high detection rates depending on uncertainties. 

with SNR 8 form lAbadie et alj ilOlOl) divided by 2.26, the correction factor for sky 
location and orientation of sources . The effect of reds hift is not included. 

^ _Dh with SNR 10 form Fig. 16 iniBake r et al.l i20Q7h . 
with SNR 8 form Fig. 3 in iReisswig et al.l <200^ . 


4i?sch = 4 • 2G'(mi + m2)/c^ because the PN approximation is 
not valid in this regime any more. 

For instance, Fig. [B] shows the time evolution of semi-major 
axis and eccentricity for a \QMq BH-BH binary with PN approx¬ 
imation. The initial semi major axis and eccentricity are 10~^ AU 
and 0.9, respectively. The solid line represents the time integration 
of Peters formula, equations 08t and ( I39t . There is a good agree¬ 
ment between the results of simulation with only 2.5PN term (open 
circle) and integration of Peters formula. On the other hand, the 
merging time of simulation with all terms up to 2.5PN (filled cir- 
cl e) is significantly small er than that of Peters formula as reported 
in iBerentzen et al.l <20091) . In case of this binary, it takes less then a 
half hour for merging (i.e., ri2 ^ 47?sch). 

In our simulations, most of pairs of close encounters have not 
been perturbed by other nearby stars. Thus, it is possible to sepa¬ 
rate the two-body motion with PN correction from the main loop 
of simulations. We, here, use a TOvQcode only for a KS two-body 
motion written by S. J. Aarseth for convenience of exploration of 
the evolution of merging binaries. The PN implementation men¬ 
tioned earlier is also adopted in the TOY code. For given semi-major 
axis and eccentricity, we simulate the orbital evolution of binaries. 
Fig.U^shows the merging time of typical BH-BH binaries in equa¬ 
tions GUI and lli with different velocity dispersion of systems. 
The merging time is estimated by two-body simulations with all 
PN corrections. A star symbol is showing the result of a galactic 
nucleus in Milky-Way-like galaxy. In all cases, the merging times 
are less than a year. Binaries formed by GR capture, therefore, will 
merge immediately. 

The GW wavefor ms of coalescing binaries have been studied 
by many authors fe.g.. lLincoln & Wilil99d : lKidderlll995h . As the 
perturbation of flat-space metric, in ca n be expresse d by (for 
more details, see equations 3.21 and 3.22 in lKiddenll995l) 


= 


2/i 

"d 


Q^J p0.5g*, ^ p gy^^ 


+ +Q%) + P^Q% + ... 


(41) 


where /r is the reduced mass, D is the distance from the source 
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Figure 16. Merging time of typical IOMq (see Olfor details) BH-BH 
binaries with different velocity dispersion of embedded star clusters. The 
merging times are obtained from two-body simulations with all corrections 
up to 2.5PN. The range of velocity dispersion is 50 to 400 km/s, which 
is correspond to the range of the mass of SMBH from 5 x IO^Mcd to 
2x10^ Mq according to the Mmbh — cr* relation from iTremaine et 51 
<2002l) . Star symbol represents the binary merging time in the Milky-Way¬ 
like galaxies. In all cases, the merging times are smaller than a year. 


to the detector, Q’’^ is the time derivative of quadrupole moment 
tensor, P^ is the PN corrections with order of n, and SO, SS and 
TT denote spin-orbit coupling, spin-spin coupling and transverse- 
traceless gauge, respectively. Here G = c = 1 is used. Since we 
are interested in the aspects of GW rather than the exact waveforms, 
we take the leading order of 




D 


i i i 7 

V V - n n-^ 

r 


(42) 


with 


=2 


i 7 i 7 

V V - n 

r 


1 


http://www.ast.cam.ac.uk/~sverre/web/pages/nbody.htm 


(43) 
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X (AU) 

Figure 17. Relative orbital motion of a BH-BH binary with the initial semi¬ 
major axis 0.153 AU and eccentricity 0.99989, as a representative of typical 
BH-BH binaries in Milky-Way-like galaxies. The orbits are very eccentric 
at the beginning. The perihelion is shifted counterclockwise due to the IPN 
and 2PN terms, and the orbit shrinks with time due to the GW emission. 

where t;* and n* are the relative velocity and the normal vector of 
the relative position, respectively. It is well known that GWs have 
two polarization -|- and x and waveforms are the linear combina¬ 
tion of these two polarizations. If we assume that the orbital plane 
lies on the xy plane initially in the source coordinate, and the angle 
between the direction to the detector N and the an gular moment um 
J is 0, the polarizations h+ and hx are given by l lKiddenll995h 

= i ^ cos^ -b sin^ — sin , (44) 

hx = cosOh^^ - sine/l^^ (45) 

Now we provide a waveform of a typical lOM© BH-BH bi¬ 
nary coalescence in a Milky-Way-like galaxy for an example. The 
semi-major axis and eccentricity after GR capture are 0.153 AU 
and 0.99989 from the equations 01 and ( I20l (. respectively. In Fig. 
[m the relative motion of BHs on xy plane is shown. Due to the 
IPN and 2PN terms, the position of perihelion is shifted counter¬ 
clockwise. In addition, by emitting GWs, the orbit shrinks more 
and more with time. Fig. [TS] shows the waveforms for this BH-BH 
binary coalescence. For simplicity, we assume that the axis of an¬ 
gular momentum is aligned with the direction to the detector (i.e., 
face-on view, 0 = 0). In Fig. IlSl a). the waveform of + polar¬ 
ization during whole evolution is presented. The merging time is 
about 8 days. Interestingly, the waveform is burst-like at the begin¬ 
ning, and it takes more than 2 days for the first burst after capture. 
The detailed waveforms h+ and hx at this moment are shown in 
lisl e) an d (d). These waveforms are similar to those of eccentric 
orbits in lAbramovici et ^ l ll992l) . Fig. [Mb), (e) and (f) show the 
waveform in the last minute, the detailed view of h+ and hx a half 
minute before merging, respectively. In this stage, the orbit is much 
circularized compared to the beginning, and the orbital frequency is 
about 10 Hz. At the moment of coalescence, the orbital frequency 
becomes few hundreds Hz which is the detectable frequency by 
ground-base GW detectors. 


One of the important questions with regards to the GR cap¬ 
tured binaries is the typical value of eccentricity when they enter 
the detectable band of the ground-based detectors. We computed 
the eccentricity when GW frequency becomes 30 Hz as a function 
of initial orbital parameters as shown in Fig. [M using the post- 
Newtonian approximation with 2.5PN term. The diagonal line is 
the locus of the typical (ao, eo) plane for different velocity disper¬ 
sion, and we computed the eccentricities at 30 Hz for up to 10% of 
peak in probability distribution function along the initial eccentric¬ 
ity for a given velocity dispersion. Different colors represent differ¬ 
ent eccentricity at 30 Hz. As can be seen from this figure, the GR 
captured binaries from typical NC with a* « 100 km/sec become 
nearly circular when they enter the LIGO/Virgo band. Only those 
from very massive NC with high velocity dispersion are likely to 
remain eccentric up to the high frequency band. Note that we have 
computed the orbital evolution only up to the pericentre distance of 
4 Rsch and some extreme binaries will not emit GW signal higher 
than 30 Hz. This could be an artifact of using post-Newtonian ap¬ 
proximation. Eventually the binaries in this regime should emit 
gravita tional waves at h igher frequencies than 30 Hz. We simply 
refer to lEast et"^ ( 1201 3h for the derivation of more realistic deriva¬ 
tion of the waveforms for eccentric binaries and their detectability 
with the ground-based GW detectors. 


8 SUMMARY AND CONCLUSION 

We have generated A^-body realizations for nuclear star clusters 
(NCs) located at the centre of galactic bulges and hosting a mas¬ 
sive black hole (MBH). In our simulations, the surrounding bulge 
is considered as the external potential well which makes the ve¬ 
locity dispersion of the embedded star clu ster isothermal si nce a 
deep potential well behaves like a heat bath dYoon et al.l201lh . The 
MBH, in the same manner, is also modeled as a point-mass poten¬ 
tial but growing with time to ensure the adiabatic adjustment of the 
stellar system. Consequently, our A^-body realizatio ns have a stel¬ 
lar density cusp (p ~ r~^ ^^t fSahcall & Wollll 19761) and Keplerian 
velocity dispersion within the radius of influence. In addition, the 
overall velocity structure is similar to observations of the star clus¬ 
ter at the centre of Milky Way (e.g.. lSchodel et alj|2009h . Strictly 
speaking, however, these star clusters are not in equilibrium but 
expand continuously since the MBH can generate kinetic energies 
by the interaction with stars in the cusp. The slopes of density and 
velocity dispersion profiles of our A^-body realizations are slightly 
shallower than those of theoretical expectations. 

This environment of NCs is a good laboratory for gravita¬ 
tional wave (GW) sources. In order to investigate GW event rates 
in NCs, we have collected the orbital information of close encoun¬ 
ters (i.e., semi-major axis and eccentricity) in our A^-body simu¬ 
lations. While most of binaries are disrupted by the strong tidal 
field from MBH, there can be many hyperbolic encounters of stel¬ 
lar mass black holes (BHs) whose pericentre distances are suffi¬ 
ciently small to radiate GWs efficiently due to the high density and 
velocity dispersion at the vicinity of the MB H and the high number 
fraction of BHs due to the mass segregation dHopman & Alexandej 
I 2 OO 6 I : l^’Learv et Ml2009h . When the energy loss by gravitational 
radiation (GR) is greater than the orbital energy, two BHs make a 
binary and merge quickly because of the small separation and large 
eccentricity after capture. Thus, the capture event rate corresponds 
to the merger rate. The capture happens most frequently near the 
half mass radius rather than within the radius of influence. Thus, 
our investigation of GR capture event rates is still valid although 
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Figure 18. Waveform of BH-BH binary coalescence for the same binary in Fig ll7l (a) for whole stage. The merging time is ~8 days. The waveform is like 
a burst, initially, (b) hy for las t minute. The waveform is much sinusoidal at this time. (c,d) /ry and h x for the first burst as marked in (a). They are similar 
to those of eccentric binaries in lAbramovici et ^il992h . (e,f) /ly and /ix at a half minute before coalescence as marked in (b). At this time, the frequency is 
~ 10 Hz. 
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Figure 19. Eccentricity of BH binaries at 30 Hz which is the typical low 
frequency limit of ground based detectors, in the initial (ao, eo) plane. The 
eccentricity at 30 Hz is expressed by the color and the diagonal solid line 
represents the typical values of (aoiGo) for different cr*. We limited the 
range of calculation down to 10% of the probability peak of the initial prob¬ 
ability distribution function of the initial eccentricity for a given velocity 
dispersion marked at the central solid line. Also we show the locus of 50 % 
of the peak as broken lines. 


our models can not precisely realize the cluster inside of the ra¬ 
dius of influence. By counting the number of GR capture events, 
we have built scaling relations of merger rates for a NC as a func¬ 
tion of the mass of the MBH and the velocity dispersion of the star 
cluster. As the result, the merger rate for a Milky-Way-like galaxy 
is ~ 10“^°yr“^ proportional to the mass ratio of MBH to the star 
cluster. 

From the Mmbh — cr* relation le.g.. iTremaine et^l2002h . 
the merger rate becomes a function of the m ass of MBH only. 
By us ing realistic mass function of MBHs (e.g.. lAller & Richstonel 
l2002h . we have determined the merger rate density per unit vol¬ 
ume. Then, the detection rates can be expressed with the merger 
rate density and the size of cosmological volume covered by GW 
detectors if we simply assume the uniformness of merger events 
over cosmic time and volume. We have obtained the expected de¬ 
tection rates 0.02 ~ 14yr“^ for advanced LIGO d epending on 
the maximum horizon distances from different studies jBaker et al.l 
l2007t l^isswig et al.ll2009l : I Abadie et alj|201^ and the mass ratio 
of MBH to the star cluster 0. 1 and 0.2. This estimate can cover 
those of lO’Learv et al.l ( l2009l) who suggested the detection rate 
of ~ 5 — 20yr“'^ of 5 to 15 Mq BH-BH binary coalescences 
in galactic nuclei for advanced LIGO. However, this study still 
have some limitations; (1) It is necessary to consider realistic mass 
function for BHs instead of assuming the mean mass fraction of 
BHs. (2) There is a relation b etween the mass of MBHs and the 
NCs jGraham & Spitleill2009h . However, we fixed the mass ratio 











































16 J. Hong, & H. M. Lee 


of MBH to the star cluster to 0.1 and 0.2. (3) In order to determine 
the maximum horizon distance for BH-BH binary mergers, the pre¬ 
cise signal-to-noise ratio calculation is needed. 

We have investigated the statistics of coalescing BH-BH bi¬ 
naries and found that the typical semi-major axis and eccentric¬ 
ity of these binaries are related to the velocity dispersion of the 
system. We also have implemented the post-Newtonian (PN) ap¬ 
proximation on the two body motions up to 2.5PN orders. With 
a given set of semi-major axis and eccentricity, we calculated the 
two-body motion under the PN approximation and the waveform 
of GW emission. The merging time is about a few hours for a 
typical BH-BH binaries in a Milky-Way-like galaxy. We also es¬ 
timated the eccentricities of GR captured binaries when they enter 
the LIGO/Virgo band (assumed to be 30Hz). Most of the binaries 
formed in relatively low velocity environment (cr, ~ 100 km/s) 
would become nearly circular while those formed in very high ve¬ 
locity environment are likely to keep significant eccentricity when 
they cross the 30 Hz. Considering the fact that the low velocity 
galaxies contribute more to the total merger rate (cf. equation [35]l, 
the GR captured binaries are likely to be indistinguishable from 
those of different origin in the sense the eccentricity effect is nearly 
negligible. 
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APPENDIX A. POST NEWTONIAN EQUATION OE 
MOTION IN CENTER OE MASS ERAME 

Here, we present the post Newtonian (PN) equation of motion of 
binary sy stem in the centre of mass frame up to 2.5P N order by 
follow ing iMora & Willi l l2004ll (also see Appendix in iBrem et S] 
|2014) . Because we assume non-spinning BHs, we do not consider 
spin terms in this paper. For the beginning, we borrow notations 
from iMora & Willi ( l2004ll as 


m = mi + m 2 , 

V = V2 - Vl, 

r = r2 - ri, 
n = r/r, 

r; = ( 77111712 ) / {mi + m 2 )^ (46) 

where mi and m2 are the mass of stars, r 1, r2, vi and V2 are the 3- 
dimensional positions and velocities, and 77 is the symmetric mass 
ratio. As mentioned in the text, the PN acceleration can be consid¬ 
ered as a perturbation and added to the gravitational acceleration as 


SL = &„ + Bipn = —^n + ^{An + Bv) (47) 

where A and B are the PN coefficients related to the relative posi¬ 
tion and velocity, respectively, where 


Ai = 2(2 + -q)— - (1 -f ‘iq)v‘^ -|- , 

r A 


A2 = - j(12 + 29??)^ - q{Z - 4q)v^ - ^77(1 

+ ^ 7(13 — 4q)—v^ -F (2 -F 2577 -F 2q^) — r^ 

2 r r 

+ ^7(3 - 477)77^7^^, 


, 8 m . f n m 

A 5/2 = TV—r( - - 

5 r V 3 r 


and 

Bi = 2(2 - q)r, 



3q)r‘^ 


(49) 

(50) 

(51) 


B 2 = —]-{4: + 41q + 8q^)—r+l;q{l5 + 4q)v^r—^q{3 + 2q)r^, 
2 r 2 2 

(52) 


8 m 




(53) 


where r is the first time derivative of radial distance defined as 
r = r ■ v/r. Then the coefficients A and B are given by the sum¬ 
mations of coefficients for different PN order divided by the speed 
of light c, Ailf?''. No te that the sign of B is opposite of that in 
[Blanchet & Iveil ( l2003h . 

Since the NBODY code uses 4th-order Hermite integrator, we 
have to have the first time derivatives of accelerations as similar to 
the accelerations 


a — fin + fipn. (54) 

The derivative of PN acceleration fipn can be expressed as 

TTl TTl 

fipn = —2— 7 ^(An+Bv)-F;^(An-FBv-FA(v/r —nf-/r)-FBa) 

(55) 

where A and B are the time derivatives of the coefficients A and 
B, which are the summations of 


Ai = —2(2 -F Tfj—r — 2(1 -F 377)v ■ a -F 3qrr, 

3 

A 2 = - (12 + 297]) ^A!—r — 4r;(3 — 47])v^v ■ a 
2 

- y 7 (l - 3q)r^r + ^qil3 - 477 )^ ^ 2 v ■ a - 

-F (2 -F 2577 -F 2q^)'!^(^2rr - 
-F 377(3 — 477 )(v • ar^ + v^rr), 


A5/2 = iq^[r--][^^+3v 


8 m . 

+ ■^n—r\ - 
5 r 


n m . 
3 r‘‘ 


and 

Bi = 2(2 - q)r, 

B 2 = -1(4 -F 4177 + 87 ")^ - A) 

-F ^^(15 -F 4?7)(2v ■ ar -F v^f) — ^t7(3 -F 2q)r^r, 


(56) 


(57) 


(58) 

(59) 


(48) 


(60) 
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B 


8 m . 

5/2 = 




8 m 
5 r 


-3^f + 2v-a 


(61) 


where r is the second time derivative of the radial distance given 

by 

f = (v^ + r • a — r'^)lr. (62) 
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